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Abstract 

Computation of frequency responses for large or- 
der systems described by time- invariant state space 
systems often provides a bottleneck in control sys- 
tem analysis. In this talk we show that banding the 
A-matrix in the state space model can effectivly re- 
duce the computation time for such systems while 
maintaining reliability in the results produced. 

Introduction to the Problem 

Consider the following realization of some transfer 
function G{ju>): 

x(t) — Ax(t) + Bu(t ) , y(t) — Cx(t) + Du(t) 

where x £ IR\ u(t) £ IR m , and y(t) £ IRC The re- 
lationship of the realization to the transfer function 
is given by 

G{juj) = C(jul- A)~'B + D 

In control system design the computation of the 
frequency response plays an important role in 
frequency-based design methods. For medium 
sized problems the order of x(t) may be in the hun- 
dreds while G(joJk) must be computed for hundreds 
of values of It has not been uncommon for a 
frequency response calculation to require hours of 
CPU time. Thus, efficient and reliable algorithms 
for this computation are needed for handling large 
order systems. 

Typical Approach 

A typical approach to computing frequency re- 
sponses for state space systems is to first perform a 
state transformation on the realization to bring the 
A- matrix into some reduced form and then solve 
the appropriate system of linear equations for each 
frequency point. 


Computational Issues 

The above algorithm for computing frequency re- 
sponses involves two issues: efficiency and sensitiv- 
ity. A potential bottleneck in computing frequency 
responses is the solution of (juj^I—A) X = B for X . 
Efficient computation is accomplished by reducing 
A to some form A which allows efficient solution of 
the above equation. Another issue is that of sen- 
sitivity. The transformation process takes place in 
finite precision arithmetic and hence will change, to 
some degree, the properties of the transfer function 
which the realization represents. It is important, 
therefore, to consider the numerical properties of 
the transformation. 

Sensitivity of the Transformation 

The effect of numerical computations in the pre- 
sense of finite precision arithmetic can be treated 
in terms of sensitivity of the coefficient matrices. 
Tranformalions which do not, increase sensitivity to 
state transformations are termed well conditioned. 
Ill-conditioned transformations can and usually do 
significantly increase the sensitivity of the coeffi- 
cient matrices to small perturbations. Presense of 
this sensitivity is often an indication of a numeri- 
cally unstable algorithm. 

Efficient Solutions to (ju^I — A)X = B 

As stated, efficient solution of the above equation is 
usually accomplished by reducing A to some spe- 
cial form. Consider the case where m — 1 (i.e., 
B £ IR. 7?xl ). Then for A in general form, solution 
of the equation requires 0(n 3 ) floating point op- 
erations (or flops). For A in Ilessenberg form or 
Schur form, where A is “nearly” upper triangular, 
solution of the equation requires 0(n 2 ) flops. The 
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transformation to produce these forms is known 
to be extremely well-conditioned. For A in diago- 
nal form or Jordan canonical form, solution of the 
equation requires 0(n) flops. However, the trans- 
formation producing these forms may be very ill- 
conditioned, leading to a reduced form which no 
longer accurately represents the tranfer function 
of interest. A good compromise is to seek some 
compromise between the upper-triangular and di- 
agonal forms. This leads reduced forms in which 
A is upper-triangular and block diagonal or upper- 
triangular and banded. We feel the banded form is 
a better strategy since it is a simpler structure to 
work with. The banded matrix is characterized by 
its order and bandwidth; the block diagonal form 
is characterized by order, block sizes, and block or- 
der. Due to the simpler structure of the banded 
matrix, algorithms based on banded matrices can 
be adopted to vector hardware architectures in a 
nicer way. 

A New Banding Algorithm 

The new banding algorithm uses several steps. 
First the matrix A is reduced to real Schur, or 
quasi-upper-triangular, form A s . Then an order- 
ing algorithm is applied to order the eigenvalues 
appearing on the diagonal of A<i in a way that will 
aid the next step in producing a small bandwidth. 
The transformations associated with the first two 
processes is very well conditioned. The third step 
involves examining the properties of the eigenval- 
ues to determine a “good” bandwidth a priori . A 
“good” bandwidth is one for which the condition 
number of T is small. Next the matrix is reduced 
to banded form, A&, using a series of operations 
to eliminate off diagonal elements. The operations 
are accumulated in a matrix T. If T is found to 
be ill-conditioned, then the tolerance for Step 3 
is tighned and Steps 3 and 4 are repeated. Fi- 
nally, the matrix A& is brought to complex, upper- 
triangular, banded form using a series of Givens 
transformations. We note that the transformations 
used in Step 4 are scaled to provide reduction in 
their condition numbers. 

An Illustration 

The figure shown illustrates the banding algorithm. 
The first operation shows the effect of bringing the 


system to Schur form. After the matrix has been 
brought to Sclmr form, the matrix is analyzed to 
determine a “good” bandwidth. Here we choose 
a bandwidth of 2. The second set of operations 
shows how the algorithm reduces a diagonal of the 
matrix. The third set of operations shows how the 
remaining diagonals are eliminated to produce the 
final upper-triangular, banded matrix. 

Test Case 

The algorithm described has been coded into For- 
tran and installed into our Pro-Matlab implemen- 
tation using the Pro-Matlab MEX facility. We 
chose as a test set a set of single input, single 
output systems with state order ranging from 20 
to 80. Matrix coefficients were generated from a 
random number generator. For each case we com- 
puted 200 frequency points. The table shows times 
for the Pro-Matlab bode function versus times for 
our bodeq function. As one can see, the new al- 
gorithm reduced the computation time from 75 to 
88 percent. 

Extensions and Future Work 

The algorithm has also been applied to time simu- 
lation of linear, time-invariant systems. The band- 
ing strategy and algorithm could be extended to 
generalized state space systems. In this case, we 
would band the A and E matrices simultaneously. 
Another possible area of future work would be pro- 
duction of better banding algorthms. The current 
algorithm is bases on solution of Sylvester equa- 
tions and has a limitation: the algorithm cannot 
band systems well when all eigenvalues are very 
closely spaced. It should be possible to band these 
matrices using different algorithms. 

Summary 

In summary, we have developed a new algorithm 
for computing frequency responses of state space 
models. In this development, we have taken into 
account the two prime computational issues: effi- 
ciency and sensitivity. We showed that the algo- 
rithms worked on a test problem and was able to 
reduce computational time considerable without a 
notable cost in accuracy. Finally, we proposed that 
the banding strategy may provide further applica- 
tion in control system design. 
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Introduction to the Problem 

Consider a transfer function G(ju>) with state space realization given by 
x(t) = Ax(t ) + Bu(t), y(t) = Cx(t) + Du(t ) 
where u(t ) G IR m , x(t) € IR”, and y(t) G 1R P . 

G(ju > ) is associated with {A,B,C, D} through 

G(ju) = C(ju>I -A)~ l B + D 

Problem: desire for many (hundreds) of values of juif, 

where n < 200 (medium order systems). 
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Typical Approach 

1 Transform the system realization: 

{A,BX\D} X {A,B,C\D} := {T~ l AT, T~ l B, CT, D} 

2 For each do 

a) solve — A)X = B for X 

b) compute G{ju>f.) = CX + D 

Compuatational Issues 

• Solution of (jto k I - A)X = B must be efficient. 

This is usually the limiting factor. 

• Given the presence of finite precision arithmetic 

{A,BX,D} must be an accurate realization of G(jui). 
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Sensitivity of the Transformation 


Consider 


{A+AA,B+AB,C+AC,D+AD} 5U {A+AA,B + AB,C+AC,D + AD} 

Then we have 


k(T) 
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where k(T) := ||T||||T- 1 || and || • || is some consistent norm. 
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Efficient Solutions to (jw*./ - A) X = B 


• A a general matrix 0(n flops/w^ 

• A an upper Hessenberg or Schur matrix => 0(n 2 ) flops/u;j. 

• A a matrix in diagonal form => 0(n) flops/u>£ 

• A a matrix in Jordan canonical form => (9(n) flops/a;*. 

• A a block-diagonal matrix =>• O(n^) flops/u.^., where 1 < / < 2 

• A a banded matrix => 0(n • 6u>) flops/w^ where bw is the bandwidth, 0 < bw < n 


75 



JPL 


A New Banding Algorithm 

Given a full general A matrix, produce A 5 where A(, is upper triangular and banded. 
The algorithm is essentially 

1 Reduce .4 to upper real Schur form, 

2 Order the eigenvalues appearing on diagonal blocks of 4 . s to produce A 0 

3 Analyze A 0 to determine a "good” bandwidth (to make k(T ) small) 

4 Uses an algorithm based on solving Sylvester equations to band A 0 , producing A f, 

5 Convert the quasi-upper triangular matrix A & to complex, upper triangular form. 


Step 4 uses transformations of the form 
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An Illustration 
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Test Case 


Single-Input, Single-Output Model with Random Coefficients 
Pro-Matlab bode function versus banded bodeq function 


pts. 

n 

bode 

(sec) 

bodeq 

(sec) 

bw 

k(T) 

200 

20 

39.9 

4.8 

3 

250 

200 

30 

81.1 

10.6 

3 

350 

200 

40 

139.8 

24.0 

13 

45 

200 

50 

211.7 

41.6 

17 

38 

200 

60 

300.6 

64.8 

20 

66 

200 

70 

407.7 

104.6 

31 

27 

200 

80 

527.1 

135.1 

32 

67 


Reduction in time from 75% to 88% 
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Extensions and Future Work 


• Application to time simulation of {.4, D , C, D } 

• Application to extended (or generalized) models: 

Ex = Ax + Bu , y = Cx + Du 

• Better banding algorithms: Banding strategy has good potential but current algorithm 
has some limitations. 
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Summary 


• Developed new algorithm for computing frequency responses of state space systems. 

• The algorithm provides a method for trading off the two computational issues at hand: 
sensitivity and efficiency. 

• The algorithm was shown to provide large saving in computational time on a set of test 
problems. 

• The strategy has some potential for other applications on medium order models. 
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